{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# II- Demo. Optimal Interpolation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "An example of simulated SSH data access is provided in the \"example_data_access_meom.ipynb\" notebook. Here, an example of a mapping technique based on a simple optimal interpolation is proposed. The notebook is structured as follow: \n",
    "\n",
    "    1) set optimal interpolation parameters,\n",
    "    2) reading of pseudo-observations,\n",
    "    3) perform optimal interpolation and,\n",
    "    4) save the results (reconstructed SSH field)\n",
    "\n",
    "\n",
    "Here, we assume a vector of observations, noted $y$ defined as:\n",
    "\n",
    "$$y = H x + \\epsilon $$\n",
    "\n",
    "where $H$ is a linear observation operator between the reconstruction grid space and the observation space\n",
    ", $x$ is the state to estimate and $\\epsilon$ is an independent observation error.\n",
    "\n",
    "The optimal interpolation consists in estimating an analysed state $x_{a}$ in combining the available observations to approximate the real state $x$:\n",
    "\n",
    "$$x_{a} = K y $$\n",
    "where $K$ is the weigth matrix defined as:\n",
    "\n",
    "$$ K = BH^T(HBH^T + R)^{-1} $$\n",
    "\n",
    "$B$ is the covariance matrix of $x$, and $R$ the covariance matrix of the error vector $\\epsilon$ ($^T$ is the transpose operator)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import xarray as xr\n",
    "import numpy\n",
    "import pyinterp\n",
    "import dask\n",
    "import warnings\n",
    "import logging\n",
    "import sys\n",
    "import os\n",
    "warnings.filterwarnings('ignore')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "logger = logging.getLogger()\n",
    "logger.setLevel(logging.INFO)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cluster = dask.distributed.LocalCluster()\n",
    "client = dask.distributed.Client(cluster)\n",
    "client"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sys.path.append('..')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from src.mod_oi import *\n",
    "from src.mod_inout import *\n",
    "from src.mod_regrid import *\n",
    "from src.mod_eval import *\n",
    "from src.mod_plot import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1) set optimal interpolation parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# OI Grid\n",
    "lon_min = 295.                                           # domain min longitude\n",
    "lon_max = 305.                                           # domain max longitude\n",
    "lat_min = 33.                                            # domain min latitude\n",
    "lat_max = 43.                                            # domain max latitude\n",
    "time_min = numpy.datetime64('2012-10-22')                # domain min time\n",
    "time_max = numpy.datetime64('2012-12-02')                # domain max time\n",
    "dx = 0.2                                                 # zonal grid spatial step (in degree)\n",
    "dy = 0.2                                                 # meridional grid spatial step (in degree)\n",
    "dt = numpy.timedelta64(1, 'D')                           # temporal grid step\n",
    "\n",
    "simu_start_date = '2012-10-01T00:00:00'                  # Nature run initial date\n",
    "\n",
    "glon = numpy.arange(lon_min, lon_max + dx, dx)           # output OI longitude grid\n",
    "glat = numpy.arange(lat_min, lat_max + dy, dy)           # output OI latitude grid\n",
    "gtime = numpy.arange(time_min, time_max + dt, dt)        # output OI time grid\n",
    "\n",
    "# OI parameters\n",
    "Lx = 1.                                                  # Zonal decorrelation scale (in degree)\n",
    "Ly = 1.                                                  # Meridional decorrelation scale (in degree)\n",
    "Lt = 7.                                                  # Temporal decorrelation scale (in days)\n",
    "noise = 0.05                                             # Noise level (5%)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2) reading of pseudo-observations + define output folder"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "inputs = ['../dc_obs/2020a_SSH_mapping_NATL60_jason1.nc'] "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load obs dataset\n",
    "ds = xr.open_dataset(inputs[0])\n",
    "ds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define outputs\n",
    "output_directory = '../results/'\n",
    "if not os.path.exists(output_directory):\n",
    "    os.mkdir(output_directory)  \n",
    "output_oi = f'../results/ssh_reconstruction_{time_min}-{time_max}_jason1.nc'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3) perform optimal interpolation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "# set OI param & grid\n",
    "ds_oi1_param = oi_param(Lx, Ly, Lt, noise)\n",
    "ds_oi1_grid = oi_grid(glon, glat, gtime, simu_start_date)\n",
    "# Read input obs + discard a bit...\n",
    "coarsening = {'time': 5}\n",
    "ds_oi1_obs = read_obs(inputs, ds_oi1_grid, ds_oi1_param, simu_start_date, coarsening)\n",
    "# Run OI\n",
    "for it in range(len(gtime)):\n",
    "    oi_core(it, ds_oi1_grid, ds_oi1_param, ds_oi1_obs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4) save the results (reconstructed SSH field)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ds_oi1_grid.to_netcdf(output_oi)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
